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Abstract 



We present a time-frequency method to detect gravitational wave signals in inter- 
ferometric data. This robust method can detect signals from poorly modeled and 
unmodeled sources. We evaluate the method on simulated data containing noise 
and signal components. The noise component approximates initial LIGO interfer- 
ometer noise. The signal components have the time and frequency characteristics 
postulated by Flanagan and Hughes for binary black hole coalescence. The signals 
correspond to binaries with total masses between 45M Q to 70M Q and with (optimal 
filter) signal-to-noise ratios of 7 to 12. The method is implementable in real time, 
and achieves a coincident false alarm rate for two detectors 1 per 475 years. At 
this false alarm rate, the single detector false dismissal rate for our signal model is 
as low as 5.3% at an SNR of 10. We expect to obtain similar or better detection 
rates with this method for any signal of similar power that satisfies certain adiabatic- 
ity criteria. Because optimal filtering requires knowledge of the signal waveform to 
high precision, we argue that this method is likely to detect signals that are unde- 
tectable by optimal filtering, which is at present the best developed detection method 
for transient sources of gravitational waves. 
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I. INTRODUCTION 



According to Thorne [[[]], gravitational wave (GW) astronomy "will create a revolution in our 
view of the universe comparable to or greater than that which resulted from the discovery of 
radio waves." He further asserts that "when gravitational waves are finally seen, they will come 
predominantly from sources we have not thought of or we have underestimated." It follows that 
GW data analysis tools should include detection methods for poorly modeled and unmodeled 
signal waveforms. 

However, at present the only well-characterized method being widely implemented for the de- 
tection of GWs from burst or transient sources is Wiener's optimal filter [j2p. This is a natural 
choice for sources whose signal waveforms are theoretically well modeled, because it is the opti- 
mal linear detection algorithm for such waveforms [Bp. Unfortunately, the effectiveness of optimal 
filtering is greatly reduced by errors in the predicted signal waveform. Furthermore, even small 
errors in GW source modeling can lead to large cumulative errors in the predicted waveform [Q] . 
Optimal filtering is therefore a poor technique for inadequately modeled (or unmodeled) sources. 
In fact, only two potential GW sources, binary inspiral and black hole quasi-normal ringdown, are 
thought to be adequately modeled for this method to work. Clearly, a method whose effectiveness 
is only weakly dependent on (or perhaps independent of) prior knowledge of the signal is needed. 
We call such methods "robust." 

One class of robust techniques widely used for signal analysis is time-frequency (TF) methods 
(cf. [£]]). The central idea is straightforward: one simultaneously decomposes the data in two 
bases, time and frequency. The resulting distribution, pit, f), represents the energy of the data 
waveform at time t and frequency /. 

Time-frequency methods are well suited to analyzing interferometric gravitational wave data. 
Since interferometers are broad band instruments, the noise energy will be distributed throughout 
the TF plane. GWs, on the other hand, are caused by bulk motions of matter and energy, and their 
spectra tend to be peaked about characteristic frequencies determined by the source dynamics. GW 
signals may therefore be identified as ridges in the surface p(t, f). While preliminary investigation 
of time-frequency methods for GW data [|6|-|T0|] show promise, to date we know of no complete 
TF detection method that has been implemented and evaluated. 

This article describes and evaluates a time-frequency method for interferometric GW data 
analysis. The method has three steps: 

1. transform the interferometer data into a Wigner-Ville TF distribution pit, f) (described in 
Section H), 



2. search for ridges in pit, f) using Steger's algorithm (described in Section p^), 

3. threshold on length of ridge to eliminate false alarms. 

We find this algorithm can reliably detect weak signals in simulated data with minimal assumptions 
about the signal. 

Our evaluation of the TF method (Section |V|) consists of estimating false alarm and false dis- 
missal probabilities for a variety of signals. False alarm probabilities are estimated by applying the 
method to a large number of data segments containing only simulated initial LIGO detector noise 
and simply counting detections. For false dismissal probability estimates, the same procedure 
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is performed on data containing both simulated initial LIGO noise and simulated intermediate- 
mass (22.5M Q — 22.5M to 35M — 35M Q ) binary black hole coalescence (BBHC) waveforms. 
These waveforms are an appropriate testing ground for robust methods, since they are probably 
dominated by the poorly understood merger signal. Furthermore, the recently discovered possible 
evidence of "middleweight" black holes [|TT],|T^p implies that these may be important sources of 
GWs. Flanagan and Hughes [[T3|] have estimated the duration and frequency band of the merger 
signal, and our test waveforms are constructed to conform to these estimates (see Appendix |S| for 
details). Thus, while our model waveform is not the correct signal for a BBHC, we believe it has 
the correct general characteristics for our evaluation. 

The results of our evaluation are promising. We obtain a false alarm rate of one per 3.4 hours 
in a single detector, or approximately one per 475 years for coincident detection in the two LIGO 
4 km detectors. Here, coincident means that signals are detected in both detectors within a certain 
time interval, or coincidence window. The coincidence window is taken to be 0.01 seconds (the 
light travel time between the LIGO 4 km detectors), since that is the maximum time interval 
between the arrival of a single gravitational wave (which travels at the speed of light) at these 
two detectors. False dismissal rates vary with the signal-to-noise ratio (SNR) (as measured by 
an optimal filter, cf. [0]) and binary mass, but as an example we find that approximately 3% of 
signals are lost at an SNR of 1 1 for a 3OM — 3OM BBHC. Graphs showing false dismissal rates 
as a function of binary mass for a range of SNRs and as a function of SNR for a range of binary 
masses are presented in Section [V|. This TF method is also computationally efficient; we are able 
to analyze data in these simulations at about twice the acquisition rate of the simulated data. 

While we are encouraged by these results, there is much more research to be done. We have 
restricted our attention in this paper to a single TF distribution (Wigner-Ville), a single ridge de- 
tection scheme (Steger's algorithm), and a single thresholding scheme (ridge length). For each of 
these there are a variety of choices available, and these choices need to be explored and evaluated. 
We have also restricted our attention to only BBHC signal models, and only for a limited binary 
mass range. While we feel that these are priority targets for such a robust search method, there are 
other sources which should be investigated. A more complete description of these and other issues 
for future research are presented in Section W\. 



II. THE WIGNER-VILLE DISTRIBUTION 

A. Time-Frequency Distributions 

The central idea of TF methods is to convert time domain data, h(t), into a time-frequency 
distribution (TFD) p(t, f). Here, the variable / labels the basis vectors in an alternative basis 
which spans the Hilbert space to which h{t) belongs (i.e. alternative to a basis labeled by t). This 
alternative basis is usually the frequency basis associated with the Fourier representation of h(t), 
although it need not be (e.g. when wavelet bases are used). Likewise, the meaning of p(t, f) 
varies, although it is usually associated in some way with the squared magnitude of h at time t and 
frequency /. 

Methods for constructing p(t, f) can be divided into two broad categories [jTJ]]: atomic decom- 
positions such as windowed Fourier transforms and wavelet transforms and bilinear distributions 
of the form 
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p(t,f) = J j h(u)h*(v)W(u,v;tJ)dudv (2.1) 
where * denotes complex conjugation and the kernel W satisfies 

W(u,v;t,f)dtdf = 8(u-v). (2.2) 



Both types of TFD might in principle be used to detect gravitational waves. 

The choice of a suitable TFD is governed by two requirements. The first is that the TFD must 
have good time-frequency localization properties. By this we mean the following: if a signal writ- 
ten in the form h(t) = A(t) cos(p(t) (with <p(t) a continuously differentiable monotonic function) 
satisfies the adiabaticity conditions 

A(t + 5(tj)*A(t), (2.3) 
<p {t + 6{tj) (2.4) 



where S(t) is defined by <p (t + 5(t)j = <p(t) + 2ir for all t, then it should give rise to a distribution 
p(t, f) which has support only in a small neighborhood of the curve / = (p(t). In other words, 
if the signal has a well defined frequency at each time, then the TFD should reflect that. This 
is required for a ridge detection algorithm to work well with the TFD. The second requirement is 
that the TFD be relatively inexpensive to compute. This makes real-time analysis of interferometer 
data feasible. 

Goncalves, Flandrin and Chassande-Mottin have investigated the localization properties 
of various TFDs when applied to binary inspiral chirp waveforms. They conclude that bilinear 
distributions have superior localization properties, and we therefore focus our attention on those. 
On the other hand, bilinear distributions can be computationally expensive, depending on the 



choice of kernel W in fl2.ip. The question, then, is whether there is a kernel for which p can be 
calculated efficiently. 



B. The Wigner-Ville Distribution 

There is at least one choice of kernel which leads to a computationally efficient bilinear distri- 



bution with good localization properties: the Wigner-Ville distribution (WVD) [ |15|J16| ] 



p(t ' /} = £L h * {* ~ i) h v + 1) et27TfTdT - (2 - 5) 

We will use the WVD , but note that it exhibits two features that might at first glance seem unde- 
sirable. 

Both features are illustrated by considering a purely sinusoidal signal, h(t) = sm(2nf t). The 
WVD for this h{t) is 

pit, f) = i [5(f - / ) + 5(f + f )- 2 cos(47r/ f)*(/)] . (2.6) 

First, note that p(t, f) does not have support only at the sinusoid's frequency / = ±/ , but also 
has an "echo" at / = 0. The sum of two sinusoids h(t) = sin(27r/ii) + sin(27r f 2 t), has further 
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echoes at / = ±(/i ± f^)- These echos are due to the bilinearity of the distribution. Second, 
observe that this p(t, f) attains negative values at / = whenever cos(47r/ £) is positive. Thus, it 
is not entirely correct to think of the p(t, f) as the squared magnitude of the signal at a give time 
and frequency. 

Despite this behavior, the WVD is well suited to our method. We wish to look for ridges in 
the TF plane. Negative values away from these ridges are of no consequence. In fact, throughout 
this paper we will adopt the convention of setting negative values of p(t, f) to 0. Also, since we 
are interested only in detection here, and not in extracting information about the signal, it does 
not matter how many ridges there are in the TFD, it matters only that their existence be highly 
correlated with the presence of a signal and that they be sharp. The WVD satisfies both these 
criteria. 



C. Discrete Implementation 

While the discussion above has been couched in the language of functions of continuous vari- 
ables, in practice it is necessary to calculate WVDs from discretely sampled data hj = h(JAt). 
This appears to presents a minor dilemma, since (|2.5[) contains expressions of the form h(t — r/2), 



which when discretized become hj_ k / 2 - For interferometer data, this issue is easily dealt with, 
since the data will be significantly (~ 10 times) oversampled. It is thus a simple matter to dec- 
imate these data to twice the time resolution required. We resample our simulated data so that 
hj = h(2jAt) accordingly. 

A second issue (also present in the continuous case) is that in practice h(t) is known only in 
a time interval < t < T. This means that h(t — r/2) is undefined for r/2 > t, and likewise 
h(t + r/2) for r/2 > T — t. In order to calculate Q2.5| ) one must assign values to h{t) in the entire 
time interval — oo < t < oo. We therefore take h(t) = in the intervals t < and t > T. 

Having resolved these issues, it is straightforward to calculate the discrete analog of (|2~3|) 



N/2 

Pjk = ^2 h(j-e/2) h(j + e/2) e 2mM ' N , (2.7) 

l=-N/2 

where N = T/ At. This discrete distribution can be treated as a digital image, with j and k 
respectively denoting the horizontal and vertical pixel number. The value of pjt denotes the gray 
scale value, or equivalently height (i.e. z value), of the image in that pixel. 



It is easy to estimate the computational efficiency of the discrete WVD from (|2.7|). For each 
value of j, one does a single multiplication (of negligible cost) and a Fourier transform. The Fast 
Fourier transform costs N\og 2 (N) floating point operations, and must be done for each of N/2 
possible values of j. Thus, the cost of calculating the WVD is approximately (N 2 log 2 N) /2. 

Finally, note that for a data set with N samples the resulting WVD has N 2 /8 pixels (N/2 
time intervals by N/A positive frequency bins). Images from large data sets therefore quickly 
become unwieldy; for example the data sets we used contained 4096 floats, leading to an image of 
2048 x 1024 pixels. We therefore averaged over 4 pixel intervals in time and 2 pixel intervals in 
frequency to obtain a more wieldy image size of 512 x 512. 

We end this Section with an example of a discrete Wigner-Ville TFD. Figure [1] shows two 
WVDs: one of simulated initial LIGO noise and the other of a simulated BBHC waveform em- 
bedded in that same noise (see Section |V A| and Appendix [A| for details). 
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FIGURES 




time (sees) time (sees) 

a) b) 
FIG. 1. The WVDs of a) simulated initial LIGO noise and b) a simulated signal and noise (optimal 

filter SNR=8). In both cases, the data is over-whitened (see Section |v|). In this figure, the darkness of a 

pixel corresponds to the magnitude of p(t, f) in that pixel, with darker pixels having higher values. Thus, 

a dark curve corresponds to a ridge in the distribution. In b), the signal, which is a model binary black 

hole coalescence waveform, is the dark structure in the bottom center (at about 0.2 sec) of the figure. The 

additional structure in the figures is self interference due to the bilinearity of the Wigner-Ville distribution. 



III. STEGER'S RIDGE DETECTION ALGORITHM 



A gravitational wave signal in interferometer data hit) should produce a ridge in the TFD 
p(t,f). Therefore, to detect GWs, a ridge detection algorithm (or equivalently line detection 
algorithm if p(t, f) is represented as a gray-scale map as in Fig. P is required. Fortunately, there 
are a number of ridge detection algorithms in the digital image processing and computer vision 



literature from which to choose QT7Q. 

We use Steger's second-derivative hysteresis-threshold algorithm [p"8[]. The essential idea of 
this scheme is simple. A ridge in a surface will have high curvature (second derivative of p(t, /)) 
in the direction perpendicular to the ridge. Furthermore, the first derivative will vanish at the top 
of the ridge, since it is a local maximum. Thus, ridges are identified as contiguous sets of point at 
which p(t, f) has a high-curvature local maximum. 

Steger's algorithm, which identifies ridges in this way, has a number of steps. First the TFD 
p(t, f) is convolved with a 2-dimensional Gaussian smoothing kernel, 



K(tJ) 



27TCT 2 



(^) 2 + ( 



t 
<rfp 



(3.1) 



where the dimensionless scale parameter a allows the preferential detection of ridges of width <, 
at p oraf p , and the parameters t p and f p are characteristic time and frequency scales of the TFD 
(we used t p = 4At and f p = 2/T, the time and frequency resolution of the digitized p(t, /)). One 
then takes the the first and second derivatives of the convolution 
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/OO pOO 
/ p(t', f')K(t -t'J- f) dt'df, 
-oo J —oo 



(3.2) 



with respect to both t and /. The second derivatives are used to find the eigenvectors (n t ,rif) of 
the Hessian matrix, 



At each point in the t — f plane, the eigenvector corresponding to the largest eigenvalue magnitude, 
|a:], defines the line in the t — f plane along which the second derivative of H(t, f) obtains its 
extremal value. For points on a ridge, these vectors will be normal to the curve described by the 
ridge in the t — f plane. If the first derivative of H(t, f) vanishes in this direction, 



then the point may be on a ridge. We call such points "potential ridge points". 

If all potential ridge points were included in ridges, one would find ridges everywhere in the 
TFD due to noise fluctuations, whether a signal was present or not. To avoid this, one thresholds 
on the value of the second derivative. However, thresholding presents conflicting requirements. 
On one hand, if the threshold is too low, it will allow too many noise ridges to be detected. On 
the other hand, if it is too high, portions of a genuine signal ridge may be missed due to deviations 
below the threshold caused by noise. To improve the detection of signal ridges while suppressing 
noise ridges, hysteresis thresholding is used. This means that there are two thresholds on the 
second derivative of p, an upper threshold, which must be exceeded by any point at which a ridge 
can start, and a lower threshold, which must be exceeded by every point on the ridge. 

Finally, ridges are identified as contiguous sets of potential ridge points which meet the hys- 
teresis thresholding criteria. Isolated potential ridge points are not defined to be ridges. Note that 
small gaps in a ridge will be "smeared out" when the TFD is convolved with the Gaussian kernel. 
Thus, gaps of less than a few af p or at p are overlooked by this algorithm, which decreases the 
probability of noise fluctuations breaking a signal ridge into many smaller ridges. Conversely, the 
minimum ridge length detected is also a few at p or af p , eliminating many shorter noise ridges. 

There are some further issues involved with implementing this algorithm on a digital image. 
When t and / are replaced by their discrete counterparts, t, and f k , the distribution p jk = p(tj, f k ) 
can be viewed as a piecewise constant function, having value pj k in the rectangle [tj — (At/2), tj + 
(At/2)] x [f k - (A//2), f k + (A//2)], where At = tj - and Af = f k - f k _ lt and vanishing 
outside the image (i.e. when j < or k < and when j > N or k > N). This piecewise constant 
function is then convolved with the continuous Gaussian kernel. The convolution is equivalent to 
the summation 




(3.3) 




(3.4) 



H, 



j k 



(3.5) 



where the convolution mask K^_ m ^ k _ n ) is given by 

rtj-tm+At rfk-fn+Af 



(j—m)(k—n) I I 

Jtj-tm-At Jfk-fn-Af 



K(t,f) dtdf. 



(3.6) 
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Rather than taking discrete derivatives of H, the process can be made more efficient by using 



masks built from the derivatives of K. By integrating < J3T2| ) by parts, one may see that the same 



result can be derived by either process. However, while the derivatives of H need to be calculated 



at each step, the derivatives of K can be calculated just once. Thus, along with Q3.6[), convolution 
masks are also made with K replaced by its first and second derivatives with respect to t and /. 

Another issue addressed by Steger's algorithm is the positioning of potential ridge points. For 
a digital image, the ridge will be composed of pixels within which the directional derivative ( |3.4| ) 
vanishes. A one dimensional example illustrates the method used to determine whether ( |3.4| ) is 
satisfied within a given pixel. Denote the one-dimensional image function by f(x). Approximate 
f(x) at a point Xj by its second order Taylor series, 

/(•'•) « fi + ^(/i+i - /i-i) {x - x 3 ) + 27^a(/i-i + fi+i ~ 2/i) (* - -O) 2 - (3-7) 

where the coefficients are finite difference approximations of the derivatives of / at Xj and Ax = 
The derivative of ( |377| ) vanishes at x = Xj + dx where 



dx = -Ax ( K- 1 ) . (3.8) 

This is within the j th pixel if and only if Xj — Ax/ 2 < Xj + dx < Xj + Ax/2. In that case, Xj + dx 
is considered to be a potential ridge point and the j th pixel is a "potential ridge pixel." One then 
joins ridge pixels (as before) to find the ridge. The generalization to two dimensions is discussed 
inpl]. 



IV. THE DETECTION METHOD 

We have thus far described two of the three steps in our GW detection method: generation 
of WVDs from interferometric data and the search for ridges (signals) in them. If all ridges cor- 
responded to GW signals, the existence of a ridge would be an adequate detection statistic and 
these two steps would be sufficient. However, even after hysteresis thresholding (see Section p^), 
stochastic detector noise will lead to spurious ridges in the WVD. This may be exacerbated if the 
noise is non-stationary and/or non-Gaussian, in which case signal ridges and noise ridges might 
have very similar characteristics. A more powerful statistic can be devised by thresholding on 
ridge characteristics which are more strongly correlated with signal ridges than with noise ridges. 
Of course, the more specific the threshold is to the signal, the less robust the detection technique 
will be. 

Ridge length (i.e. the number of pixels in a ridge) is one characteristic which distinguishes 
signal ridges from noise ridges in stationary Gaussian noise. Since the noise is stochastic, it will 
produce a distribution of ridge lengths, with short ridges being more frequent than longer ones (e.g. 
see Fig. |3|). The ridge length of a signal is (approximately) determined by the signal's duration 
and frequency band. Thus ridge length is more strongly correlated with signal ridges than noise 
ridges, and setting a minimum ridge length threshold can further improve the detection statistic. 
Also, note that this thresholding is not so specific as to undermine the robustness of the TF method. 
For these reasons, we use ridge length as a threshold. 

Our complete TF detection method consists of the following steps: 
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1. 

2. 
3. 



construct the Wigner-Ville distribution pit, f) of the detector output h(t) as per Section 0, 
search for ridges in the TF map p(t, f) using the Steger's algorithm as per Section |D 
use a threshold on the length of the ridge as a detection criterion. 



The remainder of this article evaluates the performance of this method. 



V. SIGNAL DETECTION AND FALSE ALARM STATISTICS 

The statistical performance of our TF method is measured by two probabilities: the probability 
of finding a signal when there is none, ox false alarm probability and the probability of failing to 
find a signal when there is one, ox false dismissal probability. These probabilities depend on the 
details of the noise and the signal. 



A. Noise and Signal Models 

Our simulations were carried out with discretely sampled colored Gaussian noise. This noise 
was produced using the Numerical Recipes [ |19| ] routine gasdev( ) . We altered gasdev( ) to 
use the ran2 ( ) random number generator, since it produces much longer sequences of pseudo- 
random numbers than the default rani ( ) routine. We divided the unit- variance white Gaussian 
deviates produced by gasdev( ) into two equal sets, which were used as the real and imaginary 
components of the complex valued frequency domain noise data, iij. Finally, we colored the noise 

data by multiplying them by the square root of the initial LIGO noise curve, n,jJSh{\f\) — > ftj. 
The specific Sh(\f\) we use is that of Abramovici et al. [ |20| ] as generated by the GRASP routine 
noise_power( ) with input parameter noise_f ile="noise_ligo_init . dat " 0. 

The signal data was obtained by sampling a continuous waveform developed specifically to test 
this algorithm. This waveform simulates the signal from an intermediate-mass coalescing black 
hole binary. These are expected to be important sources whose gravitational waveforms cannot 
be calculated analytically in the initial LIGO sensitivity band. They are thus ideal candidates for 
a robust detection method. Our model is based on the predictions of Flanagan and Hughes JT3|] . 
We wish to emphasize that it is not intended to be "accurate" in the sense required for optimal 
filtering: rather, it is constructed to have time, frequency and energy characteristics consistent 
with the assumptions of [JT3J] . Thus, the performance of our TF method for these simulated signals 
should be indicative of its performance for actual coalescence signals. A detailed description of 
the model waveform and its derivation is presented in Appendix A typical waveform (for a 
3OM — 3OM binary system) is shown in Fig. [| 
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FIG. 2. A simulated binary black hole coalescence waveform for a 30M© - 30M© system, showing 
the three phases of the coalescence described by Flanagan and Hughes [13]. The frequency of the wave 
increases monotonically. Only the portion of the waveform with frequency > 40 Hz is shown. 



B. The Simulation Procedure 

We analyzed simulated data in segments of N = 4096 samples each. The assumed sam- 
pling frequency was f s = 9868.420898Hz (we chose the sampling frequency of the LIGO 40- 
meter prototype detector so that data from that detector could be easily analyzed). Each segment 
therefore lasts for T = N/ f s ~ 0.415 s, and the time interval between successive samples is 
At = 1/ f s ps 1.013 x 10~ 4 s. Samples in each segment are denoted by hj = h(JAt) where j 
takes the values between and 4095. We denote the Fourier transform of h(t) by h(f) and its 
discrete representation by hj = h(jAf), where Af — 1/T and j takes the value between (DC) 
and N/2 (the Nyquist frequency). 

We performed two types of simulations: one to determine false alarms probabilities (the prob- 
ability that a signal is detected when none is present) and one to determine false dismissal prob- 
abilities (the probability that no signal is detected when one is present). In every simulation, the 
data stream contained Gaussian noise, as described in Subsection |V A| . To determine false alarm 
probabilities, no simulated signals were added, so that hj = "hj. To determine false dismissal 
probabilities, simulated GW signals were added to the noise. 

As indicated in Appendix [A], the angle averaged signals were generated in the time domain. 
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We therefore took the Fourier transform of the average signal, (h BBn ) k = J2f=o l e 2nl: > k / N (h BBH )j. 
The signal was normalized so that 

^ 2 |^BBH> fe 2 

E 1 m =1. (5-D 

fe=0 Wh)k 

where (Sh)k = Sh(kAf) is the one-sided power spectrum of the noise. The frequency domain 
data stream is then taken to be, 

h k = n k + SNR(h BB}1 } k , (5.2) 

where SNR is the signal-to-noise ratio obtainable by matched filtering. 

Since the Wigner-Ville distribution is computed from the time series representation of the de- 
tector output, it was necessary to transform h k into the time domain. However, the power spectrum 
spans such a large dynamic range that significant numerical errors arise if one simply takes the in- 
verse Fourier transform of hk- We therefore "over-whitened" the data, hk, — > hk/ (Sh)k, to reduce 
the (time-domain) dynamic range. This procedure has the added benefit that it emphasizes the fre- 
quencies in which the detector is most sensitive while suppressing the frequencies where the noise 
is large. We then took the inverse Fourier transform of the over-whitened h k , and subsequently 
computed its WVD, pj k = p(jAt, kAf), as per Eq. (|23|). 

This procedure produces discrete TFDs, whose dimensions may be expressed in units of pix- 
els. Because of the undersampling involved in computing the discrete WVD (see Section ||) we 
computed p jk at 4096/2 = 2048 distinct offsets, giving the WVD a "width" of 2048 time pixels. 
Also, the transform is constructed only at positive frequencies up to the Nyquist frequency after 
undersampling. Thus, the WVD has a "height" of 1024 frequency pixels. The resulting TFDs had 
2048 x 1024 « 2 x 10 6 pixels. This was too large for efficient computation, so we reduced the 
map size by averaging over every 4 time pixels and every 2 frequency pixels. The final discrete 
WVDs had dimension 512 x 512 pixels, each pixel having area 0.0008s x 9.6Hz = 7.8 x 10~ 3 . 

Next, pj k was passed through Steger's line recognition algorithm. We used Steger's implemen- 
tation of this algorithm pi"|], which assumes that the pixels of the image (distribution in our case) 
are unsigned characters, taking integer values between and 255. To convert the floating point 
image pj k to its unsigned character analog (p b )j k it was necessary to rescale the data to fit in this 
range without saturating it. The scaling factor was calculated as follows: the maximum floating 
point value of pj k was found for each segment of noise data. The scaling factor was then chosen 
so that it scaled the average value of these maxima to the value 128. Once rescaled, the image was 
passed to Steger's algorithm. The algorithm parameters we used were a = 2 and second deriva- 
tive (hysteresis) thresholds of 10/pixel 2 and 3.33/pixel 2 . Ideally, these values would be chosen by 
some optimization procedure, however, for this preliminary study they were selected "by hand" 
after extensive numerical experimentation. 

For each map, the line recognition algorithm returns a list of ridges detected and their lengths. 
Our detection statistic was the length of the longest ridge in the map. Thus, a threshold value was 
chosen, and if this value was equaled or exceeded by the longest curve in a given map, a signal 
was said to have been detected in that map. For instance, if one chooses a threshold ridge length 
30 pixels, a signal is said to have been detected in any map which contains a ridge whose length is 
30 pixels or longer. 
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C. False Alarm Probability 



Our goal is an algorithm which has an acceptable false alarm rate Rj. In our case, this means 
determining a ridge length which is not equaled or exceeded in maps containing only noise more 
than once in every 1/Rf seconds. To compute this threshold, we simulated noise from an ensem- 
ble of statistically independent identical detectors. In our simulation, the ensemble consisted of 
approximately 1.7 x 10 6 detectors, which corresponds to analyzing 1.7 x 10 6 x 0.415 seconds or 
about 196 hours of simulated data. 

For each member of the ensemble, we computed a WVD. We then searched for ridges in the 
WVD using Steger's algorithm and determined the length of the longest ridge. Ridges were found 
in 106 of the 1.7 x 10 6 WVDs, or about 1 out of every 16 000 maps. In those maps in which ridges 
were found, the length of the longest ridge ranged from 7 to 68 pixels. The relative frequency with 
which the longest ridge had a length > £ pixels is shown as a function of £ in Figure [3[ 




FIG. 3. The relative frequency with which Steger's algorithm detected ridges with length > £ pixels in 
WVDs of simulated initial LIGO noise. Only the 106 maps in which ridges were detected are included in 
this graph. The unevenness of the graph is caused by statistical fluctuations due to small number statistics. 

Based on this simulation we chose a ridge length threshold of 30 pixels. This threshold pro- 
duces a false alarm probability of about 3.4 x 10~ 5 , or a false alarm rate of about one every 3.4 
hours in a single detector. More interesting is the rate at which false alarms occur simultaneously 
in two uncorrelated detectors. This coincidence rate is proportional to the coincidence window 
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r, i.e. the time interval within which the same signal can be seen by both detectors. A natural 
lower limit for r is the light crossing time between detectors, r c , since the time interval between 
the arrival of a single signal at the first detector and the last detector can be up to r c . However, 
other considerations may set a larger r for a given algorithm, leading to a higher coincident false 
alarm rate. 

For the LIGO 4 km detectors, r c is approximately 10 ms. As implemented, the TF method did 
not resolve time-of-arrival to this accuracy, because we did not distinguish when in the 400 ms 
data segment the signal arrived. The coincidence window is therefore 400 ms, which corresponds 
to a coincident false alarm rate of approximately one 400 ms/(3.4 hr) 2 , or about one per 12 years. 
However, recall that the time resolution of the WVD is less than 0.0008 s. Furthermore, recall 
that Steger's algorithm identifies the position of ridges in the image (in fact, this feature was 
implemented in our code although we did not use it). Thus, at no extra computational cost, one 
could set r = r c = 10 ms. The TF method can therefore easily achieve a coincident false alarm 
rate of 10 ms/ (3.4 hr) 2 , or about one every 475 years. 



D. False Dismissal Probabilities 

The false dismissal probability depends on the nature and strength of the signal in the data. 



We've estimated these probabilities for simulated BBHC waveforms (Section |V A[ ). We used 
signals with SNRs (as measured by optimal filters) in the range 7 to 14. To demonstrate the ro- 
bustness of this TF method we selected various coalescence waveforms corresponding to different 
binary system masses. We confined ourselves to the (total system) mass range 45M to 7OM . 
In this range the merger phase, for which the waveform cannot be calculated analytically, sweeps 
through the frequency band of maximum sensitivity for initial LIGO detectors, and dominates the 
detectable signal [p"3[|. These are therefore masses for which robust methods will be most useful. 
Figure^ shows false dismissal probabilities as a function of binary mass at SNRs of 10, 12 and 

14. 
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FIG. 4. False dismissal probability as a function of mass. The three curves correspond to three different 
values the optimal filter signal-to-noise ratio. With the parameters we have chosen, our method tends to 
work better for higher mass binaries, where the energy is more localized in the TF map. 

At an SNR of 14 the false dismissal rates are acceptable throughout the tested mass range. 
At an SNR of 10, however, more than 20% of signals are missed by our detection algorithm for 
system masses < 55Af . This is because as the mass of a system decreases, the proportion of the 
SNR which is attributable to the inspiral phase increases, as does the duration of the observable 
portion of the inspiral. Thus, the SNR of the signal is spread through a larger region of the TF 
map, leading to lower ridges and hence a loss of detectability. 

In Figure |5| we plot the false dismissal probabilities versus the SNR for signals from 
45M , 6OM and 7OM binary systems. Again, one sees in this Figure that high mass systems are 
more readily detected. False dismissal probabilities for every system mass decrease with signal- 
to-noise ratio as expected. 
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snr 

FIG. 5. False dismissal probability as a function of SNR for three values of the total mass M of the 
binary system for the model signal. Even at a modest SNR of 1 1 , the false dismissal rate is below 50% for 
systems with masses between 70M© and 45M . 



E. Computational Efficiency 

All our simulations were carried out on a 48 node Beowulf (parallel) computer comprised 300 



MHz Alpha 21 164 machines [|22[]. This computer can analyze 47 segments of data simultaneously 
(the 48th machine is used to coordinate the calculation). We found that we could apply our TF 
method to approximately 2 x 10 4 data points per second (about twice the sampling rate of the 
simulated data). 

For initial LIGO, the frequency band of maximum sensitivity is projected to be below a few 
hundred Hz [|2~0"|]. Sampling rates of <, 1024 Hz should therefore be sufficient to contain all signals 
detectable with this TF method. Thus, even modest parallel computing facilities could implement 
this method in real time for initial LIGO data. 



VI. CONCLUSION 
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A. Summary 



Optimal filtering has dominated the interferometric GW data analysis literature. This is be- 
cause when one seeks signals of nearly exactly known form in stationary noise with known spec- 
tral density, optimal filters are the most sensitive linear filters. However, it is clear that there are 
sources of gravitational radiation that are not understood well enough to predict their signal wave- 
forms. For such sources, matched filtering will not be a suitable detection method. It is therefore 
necessary to develop methods which can detect signals from poorly modeled or even unmodeled 
sources in interferometric GW data. 

We have demonstrated that time-frequency methods may provide a tool for detecting such 
sources. Our method, with the parameters discussed above, produced an initial LIGO false alarm 
rate of « (3.4 hr) _1 , which corresponds to coincident false alarm rate in two detectors of w 
(475 yr) _1 . Our simulations also show that we reliably detect the postulated BBHC waveforms for 
binaries of various masses injected into the data at a variety of optimal filter signal-to-noise ratios. 
For instance, we detect 80% of signals with SNRs ranging from « 8.8 for higher mass binaries 
(7OM ) to ~ 12.2 for lower mass binaries (45M ). These numbers are within a factor of two of 
those for optimal filtering of stellar mass binary inspirals, where for two detectors the detection 
threshold is usually taken to be an SNR of ~ 6 [|2~3|]. Furthermore, this method is implementable 
in real time with a realistic computing budget. 

B. Future Directions 

While this method shows considerable promise, there are still a number of issues which must 
be addressed in order to determine how useful TF methods will be in general: 

Utilization - One might use TF methods as a filter to identify data which should be analyzed for 
the presence of signals using other techniques (i.e. in an hierarchical search), for detecting 
the signals themselves, or as a method for characterizing detector noise. This paper only 
addresses signal detection. 

Choice of Algorithm - We have presented only one of many TF methods that could be imple- 
mented. We have used the Wigner-Ville distribution and Steger's ridge detection algorithm, 
however there are a number of other TFDs (cf. [§]) and ridge detection algorithms (e.g. 
Hough transforms [ |24]j25| ] or the curve and edge detection algorithms cited in the introduc- 
tion of [[T8]]) which might prove suitable. As a detection statistic, energy along ridges, total 
ridge length in a map, or "bunching" of ridges might provide better detection statistics than 
maximum ridge length. 

Optimization - We have made heuristic choices for operating values of the smoothing parameter 
a and the second derivative thresholds in Steger's algorithm. It would be desirable to have a 
process which, given a specific type of GW signal, optimizes these parameters. 

Implementation - It is unlikely that a single TF method is sufficiently robust to detect all types of 
unmodeled sources, which means that a set of TF methods should be used. It is also possible 
that TF methods are optimally implemented hierarchically. For example, one might have two 
detection statistics; one with an unacceptably high false alarm rate but low computational 
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cost and another with prohibitive computational cost but a low false alarm rate. By applying 
the second statistic only to those TF maps in which a signal is detected with the first statistic, 
one could achieve both computational efficiency and acceptable false alarm rates. 

Comparison with Other Techniques - We do not know how TF methods compare to other tech- 
niques. Signals that are optimal for TF detection produce regions of high density in the TFD. 
For our simulated signals, this high density region spans relatively short time scales 0.2 
s). They might therefore be equally well detected with much simpler time-series threshold- 



ing techniques such as the one described in [£6fl. We have investigated this possibility by 



comparing each datum used in our TF simulations to a threshold. We chose the threshold to 
produce the same false alarm rate as obtained with the TF method. This technique produced 
a false dismissal rate of about twice that of the TF method. However, it is unknown whether 
a more sophisticated time-series threshold or a different statistic such as the excess power 
statistic [ |27l ] would prove superior to the TF method. 

There is one further issue which pertains to all robust methods. In order to detect a signal 
of unspecified form in noise, the noise must be well characterized. We have assumed Gaussian 
noise in this paper. It is uncertain at this time to what extent this assumption is valid, however it is 
known that noise from the LIGO 40 m prototype contains a significant non-Gaussian component. 
The extent to which this component can be identified and removed (perhaps with the aid of TF 
methods themselves) is an issue that will only be fully resolved once the relevant interferometer 
data becomes available. 

Clearly, there is a much work to be done. Moreover, with interferometric detectors scheduled 
to begin taking data in approximately a years time, there is little time in which to do it. In an effort 
to expedite the resolution of the issues listed above, we have made our TF computer code available 
as part of the latest release of GRASP [@]. 

However, the results of this preliminary investigation are encouraging. With the resolution of 
these issues, TF methods promise to be useful tools for the detection of GW signals in interfero- 
metric data. 
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APPENDIX A: THE BINARY BLACK HOLE COALESCENCE WAVEFORM 

The signal waveforms we used in our false dismissal simulations are based on the description 
of binary black hole coalescence (BBHC) signals in Flanagan and Hughes (FH) [[13]]. They consist 
of three consecutive components: an inspiral component, a merger component, and a quasi-normal 
ringdown component. The division into these components is somewhat arbitrary. Roughly speak- 
ing: 
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The inspiral component describes the evolution of the binary system up to the "innermost 
stable circular orbit" (ISCO). This component is well modeled for BBH systems of ~ 1M 
], but there is some doubt whether it is known well enough for optimal filtering for 



binaries of total mass >, 20M Q [p9[|. 



The merger component describes the evolution of the system from ISCO until such time as 
the system can be described as a single perturbed Kerr black hole. No analytical descriptions 
of this component exist. 

The quasi-normal ringdown component describes the evolution of the system when it is 
well described as a single perturbed Kerr black hole. This component well modeled for all 



BBHC's [60,31] 



In this Appendix, we model each component separately. The full coalescence waveform is then 
simply a combination of the three component waveforms: 

r (h insp (t)), t < o, 

(h BBU (t)) = I A™^ e (t) cos^ mer s e (t), < t < t meige , (Al) 

K (0)> ^ ^* Emerge; 

An example of such a waveform for a 3OM - 3OM system is shown in Figure |2|. We use units 
in which c = G = 1 throughout the remainder of this Appendix. 



1. Quasi-Normal Ringdown 

The ringdown waveform is given in Jl3| ] as 

hf* _ ih f = *M 2iS f ( t p a ) e -W^*-*/r-H<Po. (A2) 
r 

Here, M and a are the mass and angular momentum per unit mass of the black hole, r is the 
distance to the hole, 25*2 { L i A a ) ^ s a spin-weighted spheroidal harmonic, / qnr is the characteristic 
frequency of the black hole normal mode, r is the time constant for the exponential decay of 
the mode, cp Q is an arbitrary phase factor, and t and (3 are angles that describe the orientation of 
the black hole's spin axis with respect to the plane of the detector arms (c.f. Fig. 9.2 of [JT]]). 
This waveform is quite straightforward except for the presence of the spin-weighted spheroidal 
harmonic, f3, a). 

This complication of the 2-S , 2( i > A a ) f ac t° r is easily circumvented by considering only the root- 
mean-squared (RMS) average of it over the source orientation angles 1 and (3. When averaged in 
this way, spin-weighted spheroidal harmonic obeys the trivial relationship 

f d Q \ 2 S!(L,P,a)\ 2 = JL (A3) 

Since we are not interested in modeling the waveform exactly (indeed, our point is to show that we 
can find even inaccurately modeled signals), let us for the sake of simplicity substitute the RMS 
value of (Any 1/2 for 2 S%(i, (3, a) in ( gg ) to obtain 
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(hDw) - WDM = 4^- e- 2 ^-^°. (A4) 

'Aur 



Note that (hf*) M) and (/^ nr ) (t)/3 ) are no? the (t, /3) angle averaged RMS values of hf 1 and /if r 
respectively, but are rather the result of averaging 2 5'| over (t, (3). 

From QA2[), the detector response to the quasi-normal ringing of a black hole is obtained by 



h(t) = F+{6, 0, if>)h+(t) + F x (9, 0, $)h x (t), (A5) 
where F + and F x are the beam pattern functions 

F + (9, 0, if>) = - (1 + cos 2 0) cos 20 cos 2^ - cos sin 20 sin 2ip, (A6) 

F x (9, 0, V) = - (1 + cos 2 0) cos 20 sin 2-0 + cos 9 sin 20 cos 2-0, (A7) 

and 9, and ip are angles which describe the detector orientation [JTJ]. One may average over 
detector orientations to obtain 



/ 1 rn rlTV rir 



1 ^h + (ty + h x (ty. (A8) 



Substituting ( |A4| ) into ( JA8D we get 

„4M 



- ^ = ^e"^. ,A9) 

Notice that A qnT (t) is no longer a waveform. The source angle averaging combines h + and h x to 
in such a way that only the (exponentially decaying) amplitude envelope remains. To recover the 
waveform, we multiply A qnr (t) by a cosine at the appropriate frequency, 



qnr\ 



(h(t)^} EE A**(t - t merge ) C0S(27T/ C » nr (t - Wge) + <P% 

_ _ J AiW_ o _ (t _ w)/T cos(27r/qnr(t _ w) + ^ nr) _ (A1Q) 



\/207rr 



Here, the angle brackets serves to remind us that there have been two averaging processes used in 
obtaining this result, and we have offset the time variable t by t mer ge (the duration of the merger 
component) to facilitate combining it with the other waveform components as in ( |A1[ ). 

Finally, we wish to express flA10| ) in terms of the BBH system mass. Using Eq. (3.17) of [ JT3| ] 
and the values a = 0.98 and A = 0.4 quoted in Eq. (3.21) of JT3|], we have 



/'"'= 1320Hz (^f 2 ) , (All) 

(A.2) 



qnr 



Thus, the only free parameters remaining in ( |A1U| ) are the system mass M, the time offset t mer ge 
and the initial phase of the ringdown waveform p qnr . The mass will remain a free parameter in our 
simulation, however t merg e and the initial phase are determined by the merger waveform. 
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2. Binary Inspiral 



For the binary inspiral component, h+ and h x can be conveniently written in the form [jl]] 

hftt) = -2 (£\ [7rM/ insp (t)] 2/3 (l + cos 2 l) cos^ insp (t), (A13) 
h}™*(t) = -4 [nMf insp (t)f 3 cost sin^ insp (t), (A14) 

where M = mi + m-i is the total mass of the BBH system, [i = (mim 2 ) jM is the reduced mass 
of the system, and f msp (t) and (p msp (t) are the binary inspiral orbital phase and frequency. These 
latter are given to first post-Newtonian order by [@] 



/ insp (t) 



— |e- 3 / 8 + ( 



743 11 * , 



2688 32 



coal 



M e5/8 + (^ + 5596 ") 93/8 } 



where 



77 = 
6 = 



JL 

M' 
g 
5M 



(^coal t) ■ 



(A15) 
(A16) 

(A17) 
(A18) 



<^coai an d t C oai are the phase and time at which coalescence occurs (i.e. when the frequency be- 
comes infinite). Note that inspiral waveforms are available to post-Newtonian order 2.5 [Q], but 
first post-Newtonian order will be sufficient for this crude model. 

As with the ringdown waveform, we perform an RMS averaging over the source angles in the 
inspiral waveform. The equivalent averaging in this case replaces (1 + cos 2 t)/2 in h+(t) and cos 1 
in h x (t) by their (4, 0) averaged RMS values, 



47T 



[1 + COS 2 L) 2 + COS 2 L 



2 

75' 



(A19) 



Thus, 



v5r / 



[vrMf nsp (t)] 2/3 cosy? insp (t), 
[vrM/ insp (t)] 2/3 sin^ insp (t), 



(A20) 
(A21) 



We also average over the detector orientations. Inserting ( |A20D and ( |A21[ ) into ( |A8[ ) and multiply- 
ing by cos (p msp (t) we have 



(/i insp (t)) = A msp (t) cos^ insp (t) = ^[7rM/ insp (t)] 2/3 cos^ insp (t). 



(A22) 



Again, we wish to express the waveform in terms of the system mass M. It therefore remains to 
fix the reduced mass /i, the coalescence time t coal , and the coalescence phase f coa i- We restrict our 
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attention to equal mass binaries, so that the reduced mass is therefore /x = M/4. We fix i coa i so that 
the binary system attains the ISCO frequency, /isco> at ti me t = 0, i.e. so that / msp (0) = /isco- 
Following [0, we use the Kidder- Will- Wiseman [0] value of / ISCO = (2OM /M) x 205Hz. 
We also fix (p coa ,\ so that the phase of the binary system vanishes at t = 0, i.e. so that y? msp (0) = 0. 
Note that this is not the same as choosing £ coa i = and (p coal = 0, since coalescence occurs after 
the ISCO. This completes the specification (apart from the free mass parameter) of the inspiral 
component of the BBHC waveform. 



3. Merger 

The remaining task is to model the merger waveform. While no analytic model for the merger 
exists, Flanagan and Hughes [jT3p make educated estimates of some of its properties. They as- 
sume that the merger signal contains only frequencies between the ISCO frequency, /isco> and 
the quasi-normal ringing frequency, / qnr . They also estimate that the energy carried by the merger 
component of the wave is approximately 3 times the energy carried by the quasi-normal ring- 
down component. Finally, they estimate the duration of the merger to be t mcrgc ~ 50M. We use 
these criteria, along with the requirement that h and d t h be continuous, to guide us in making a 
mock merger waveform. While the resulting waveform won't be that of a real BBH merger, it 
should resemble it enough to determine whether our TF method can detect BBH waveforms for 
intermediate-mass systems. 

The first step in our construction of a merger waveform estimate is to assume a form of 

(/i merge (t)) = A merge (t) cos^ merge (t). (A23) 

Continuity of the / merge and <9 t / mer g e with both the inspiral waveform that proceeds it and the 
quasi-normal ringdown waveform that follows constitutes four conditions on tp meTge (t): 

£ merge 

(* = 0)=/bco, (A24) 

'02 merge \ /^2 insp\ 

^ = 0)=(^^J(t = 0), (A25) 

(A26) 
(A27) 

where the last equation follows from the fact that the / qnr is constant. To satisfy these four con- 
ditions will require a frequency model with four free parameter. We use the simplest such model; 
the merger frequency a cubic function of time. The phase is therefore a quartic of the form 

insp \ 

<P mcrec (t) = /moo* + ^2-(°) t2 + ^T rse t 3 + vT rge t\ (A28) 




where ip™ CIgc and yj™ erge can be obtained by solving ( |A26| ) and ( |A27| ) and we have set yj merge (0) 



0. Once the merger phase polynomial has been determined, it is a simple matter to find the phase 
constant for the quasi-normal ringdown, 
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cpr = V merge (t = t merge ), (A29) 

which completes the specification of the quasi-normal ringdown component of the waveform. 
To determine the merger amplitude function, we impose similar continuity conditions 



") (t = 0)= (^ inSP ) (t = 0) ' (A31) 



A mcrgc {t = 0) = A insp {t = 0) (A30) 

/ $^merge ' 

A (t Emerge) (t Emerge); (A32) 

/Q A merge\ / Q \ 

I 7^ J ^mergc) ( ~q^A^ J ^merge); (A33) 

We again need a fitting function with at least four parameters. However, there is a further constraint 
to impose; recall that [ |T3"| ] estimates the energy of the merger to be three times the energy of the 
quasi-normal ringdown, i.e. 

/Emerge ( Q \ 2 r°° f 8 \ 2 

I ( S <*— (*)>) * = 3/ w ( S ™) *• (A34) 

This fifth condition on y4 merge (t) requires a fifth parameter, and we therefore model the merger 
amplitude with a quartic, 

A mer g e^ = ^insp^ + (0) t + A™*^ 2 + A^t 3 + A^t* , (A35) 

where the coefficients A^f 6 are determined by ( |A32p , ( |A33| ) and ( |A34[ ). This completes the 
specification of the merger waveform. 
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